The Wave Action Equation in Cartesian  (\ref{eq:bal_plane}) or spherical  (\ref{eq:bal_sphere}) coordinates is the basic
equations of the wave model. However, modified versions of these equations are
used in the model, where (a) they are solved on a variable wavenumber grid
(see below), where (b) modified versions of these equations are used to
properly describe dispersion for discretized equations in selected numerical
schemes (see \para\ref{sub:xy_prop}), and where (c) sub-grid obstacles such as
islands are considered (see \para\ref{sub:xy_prop}).


\vssub
\subsection{~Spectral discretization} \label{sec:basic_num}
\vsssub


If Eq. (\ref{eq:bal_plane}) or Eq. (\ref{eq:bal_sphere}) is solved directly, an
effective reduction of spectral resolution occurs in shallow water
\citep[see][]{tol:GAOS98b}. This loss of resolution can be avoided if the
equation is solved on a variable wavenumber grid, which implicitly
incorporates the kinematic wavenumber changes due to shoaling. Such a
wavenumber grid corresponds to a spatially and temporally invariant
grid in relative frequency \citep{tol:GAOS98b}. The corresponding local wavenumber grid can be
calculated directly from the invariant frequency grid and the dispersion
relation (\ref{eq:disp}), and hence becomes a function of the local depth
$d$. To accommodate economical calculations of $S_{nl}$ and allow a good separation of swell frequencies, a frequency discretization with 
exponentially increasing increments is adopted, so that the varying frequency resolution is proportional to the local frequency, 

%----------------------------%
% Exponential frequency grid %
%----------------------------%
% eq:sigma_grid

\begin{equation}
\sigma_{m+1} = X_\sigma \, \sigma_m \: , \label{eq:sigma_grid}
\end{equation}

\noindent
where $m$ is a discrete grid counter in $k$-space. $X_\sigma$ is defined by
the user in the input files of the program. Traditionally, in most
applications of third-generation models $X_\sigma \simeq 1.1$ is used.

The effects of a spatially varying grid will be discussed for the Cartesian
Eq. (\ref{eq:bal_plane}) only. Adaptation to the spherical grid is
trivial. Denoting the variable wavenumber grid with $\kappa$, the balance
equation becomes

%------------------------------%
% general equations kappa-grid %
%------------------------------%
% eq:bal_f_grid
% eq:sigma_dot

\begin{equation}
\frac{\p}{\p t} \frac{N}{c_g} +  
\frac{\p}{\p x} \frac{\dot{x} N}{c_g} + 
\frac{\p}{\p y} \frac{\dot{y} N}{c_g} + 
\frac{\p}{\p \kappa} \frac{\dot{\kappa} N}{c_g} + 
\frac{\p}{\p \theta} \frac{\dot{\theta} N}{c_g}  =
\frac{S}{\sigma c_g} \: , \label{eq:bal_f_grid} \end{equation} \begin{equation}
\dot{\kappa} \:\frac{\p k}{\p \kappa} =
     c_g^{-1} \frac{\p \sigma}{\p d} \left (
    \frac{\p d}{\p t} + {\bf U} \cdot \nabla_x d \: \right ) -
    {\bf k} \cdot \frac{\p {\bf U}}{\p s}
\: . \label{eq:sigma_dot}
\end{equation}

\noindent

\vssub
\subsection{~Splitting of the wave action equation} \label{sec:basic_num}
\vsssub

In \ws\, Eq. (\ref{eq:bal_f_grid}) is solved using a fractional step method. 
The first step treats the temporal variations of
the depth, and corresponding changes in the wavenumber grid. As is discussed
by \cite{tol:GAOS98b}, this step can be invoked sparsely. By splitting off
effects of (temporal) water level variations, the grid becomes invariant, and
the depth becomes quasi-steady for the remaining fractional steps. Other
fractional steps consider spatial propagation, intra-spectral propagation and
source terms. Starting with version 5.10, the source term $S$ is further split into 
non-ice $S_{no~ice}$ and ice $S_{ice}$ source term. For a single model grid, the following sequence of integration  is performed 
by the {\code W3WAVE} routine:
\begin{itemize}
 \item[1.] Update of water level 
 \item[2.] Intra-spectral part 1: integration over $\Delta t_g/2$ of  $\frac{\p}{\p t} \frac{N}{c_g} + \frac{\p}{\p \kappa} \frac{\dot{\kappa} N}{c_g} + \frac{\p}{\p \theta} \frac{\dot{\theta} N}{c_g}  =0 $ 
 \item[3.] Spatial propagation: integration over $\Delta t_g$ of $\frac{\p}{\p t} \frac{N}{c_g} +  
\frac{\p}{\p x} \frac{\dot{x} N}{c_g} + 
\frac{\p}{\p y} \frac{\dot{y} N}{c_g}  = 0$
 \item[4.] Intra-spectral part 2: integration over $\Delta t_g/2$ of  $\frac{\p}{\p t} \frac{N}{c_g} + \frac{\p}{\p \kappa} \frac{\dot{\kappa} N}{c_g} + \frac{\p}{\p \theta} \frac{\dot{\theta} N}{c_g}  =0 $ 
 \item[5.] Source term integration: integration over $\Delta t_g$ of $\frac{\p}{\p t} \frac{N}{c_g}   = \frac{S_{no~ice}}{\sigma c_g}$
 \item[6.] Ice source term integration: integration over $\Delta t_g$ of $\frac{\p}{\p t} \frac{N}{c_g}   = \frac{S_{ice}}{\sigma c_g}$ 
\end{itemize}
The succession of these 6 steps is, in the limit $\Delta t_g \rightarrow 0$, equivalent to the integration of Eq. (\ref{eq:bal_f_grid}) over a global 
time step $\Delta t_g$. 

This splitting in multiple steps allows an  efficient 
vectorization and parallelization at the same time. The time splitting furthermore
allows for the use of separate partial or dynamically adjusted time steps in
the different fractional steps of the model. \ws\ makes a distinction between
4 different time steps. \label{dt_list}

\begin{list}{xx}{\itemsep 0mm \parsep 0mm \rightmargin 5mm}

\item[1)] The `global' time step $\Delta t_g$, is the common step of all the splitted sub-integrations. In that sense, 
it is the smallest time step for which a physically meaningful solution can be obtained, because all terms in the equation have been 
integrated. As a result, this is a possible time step for evaluating model output or coupling with other models, 
and, in the case of a multi-grid system,  it is the time step at which communication 
between grids is performed. In the case of a forced -- not coupled -- model, input winds and currents are
interpolated at this global step. This time step is provided by the user in the input file of {\code ww3\_grid}, but can be reduced
within the model to reach a requested input or output time.

\item[2)] The second time step is the time step for spatial propagation. This is not used for triangular-based 
grids, for which the advection step is -- in the case of explicit schemes -- adjusted internally for each spectral component. For other grid types, the
user supplies a reference maximum propagation time step for the lowest model
frequency $\Delta t_{p,r}$, assuming no currents, and no grid motion. For the
frequency with counter $m$, the maximum time step $\Delta t_{p,m}$ is
calculated within the model as

%-------------------------%
% Time stepping equations %
%-------------------------%
% eq:dtpl

\begin{equation}
\Delta t_{p,m} = \frac{{\dot{x}}_{p,r}}{{\dot{x}}_{p,m}} \Delta t_{p,r}
\: , \label{eq:dtpl} \end{equation}

\noindent
where $\dot{x}_{p,r}$ is the maximum advection speed for the longest waves
without currents or grid motion, and $\dot{x}_{p,m}$ is the actual maximum
advection speed (including current) for frequency $m$. If the propagation time
step is smaller than the global time step, the propagation effects are
calculated with a number of successive smaller time steps. This generally
implies that several partial time steps are used for the lowest frequency, but
that the highest frequencies are propagated over the interval $\Delta t_g$
with a single calculation. The latter results in a significantly more
efficient model, particularly if higher-order accurate propagation schemes are
used. Note that $\Delta t_{p,m}$ may be defined bigger than $\Delta t_g$, and
that this has potential impact in model economy for cases with (strong)
currents.

\item[3)] The third time step is the time step for intra-spectral
propagation. For large-scale and deep-water grids this time step can generally
be taken equal to the global time step $\Delta t_g$. For shallow water grids,
smaller intra-spectral propagation time steps allow for larger effects of
refraction within the stability constraints of the scheme. Note that the order
of invoking spatial and intra-spectral propagation is alternated to enhance
numerical accuracy. If strong refraction of long period swells occur, this may
result in a notable undulation of mean wave parameters. This can be avoided by
setting this time step to an even integer fraction of $\Delta t_g$.

\item[4)] The final time step is the time step for the integration of the
source terms, which is dynamically adjusted for each separate grid point and
global time step $\Delta t_g$ (see \para~\ref{sub:source}). This results in
more accurate calculations for rapidly changing wind and wave conditions, and
a more economical integration for slowly varying conditions. In order to limit the 
calculation time, a minimum time step is defined by the user. 

\end{list}

\vspace{\baselineskip} \noindent 

The following sections deal with the separate steps in the fractional step
method, and various subjects associated with this.  The main issue are covered
in \para\ref{sub:num_depth}, which addresses treatment of temporal variations
of the water depth, \para\ref{sub:xy_prop} which addresses spatial
propagation, \para\ref{sub:spec} which addresses intra-spectral propagation,
and Sections \ref{sub:source} and \ref{sub:icesource} which address the numerical 
integration of non-icea and ice source
terms.  The other sections deal with additional
numerical approaches and techniques, covering the treatment of winds and
currents (\para\ref{sub:num_w_c}), including tides (\para\ref{sub:num_tide}), calculating space-time extremes (\para\ref{sub:space_time_ext}), treatment
of ice (\para\ref{sub:num_ice}), spectral partitioning and the corresponding
tracking of wave systems in space and time (Sections \ref{sub:num_part},
\ref{sub:num_track}), and nesting (\para\ref{sub:num_nest}).
